#include "./cell.hpp"
#include "yocto/math/fcn/zfind.hpp"

double Cell:: ComputeZeroLeak( double E )
{
    const double    zeta  = __Faraday__ * E / ( __R__ * Temperature );
    const solution &S_in  = sol;        //-- assume solution @t=0
    const solution &S_out = tmp;        //-- assume outside  @t=0
    acc.ldz();
    Leak(acc, zeta, S_in, S_out);        //-- gather leaks
    double &lambda_K = acc["K+"];      //-- modify potassium term
    lambda_K *= 1.5;
    return acc.sum_zC();
}

void Cell:: ComputeSteadyEm()
{
    sol.copy( inside(0.0)  );
    tmp.copy( outside(0.0) );
    const double E_step = 1e-3;
    double E_min = -E_step;
    double E_max =  E_step;
    while( ComputeZeroLeak(E_min) * ComputeZeroLeak(E_max) >= 0 )
    {
        E_min -= E_step;
        E_max += E_step;
    }
    zfind<double>             solve(1e-7);
    numeric<double>::function zproc( this, & Cell::ComputeZeroLeak );
    Em = solve(zproc,E_min,E_max);
    std::cerr << "Natural Em= " << Em * 1000 << " mV" << std::endl;
    
}

double Cell:: EnforceZeroLeak( double Coeff )
{
    const double zeta = __Faraday__ * Em / ( __R__ * Temperature );
    const solution &Sin  = sol;        //-- assume solution @t=0
    const solution &Sout = tmp;        //-- assume outside  @t=0
    perm_coeff["K+"] = Coeff;
    acc.ldz();
    Leak(acc, zeta, Sin, Sout);        //-- gather leaks
    double &lambda_K = acc["K+"];      //-- modify potassium term
    lambda_K *= 1.5;
    return acc.sum_zC();
    
}

void Cell:: EnforceSteadyEm( double E )
{
    Em = E;
    std::cerr << "Enforcing Em=" << E << std::endl;
    sol.copy( inside(0.0)  );
    tmp.copy( outside(0.0) );
    
    double Amin = 1;
    double Amax = 1;
    while(  EnforceZeroLeak(Amin) * EnforceZeroLeak(Amax) >= 0 )
    {
        Amin *= 0.9;
        Amax /= 0.9;
    }
    
    zfind<double>             solve(1e-7);
    numeric<double>::function zproc( this, & Cell::EnforceZeroLeak );
    perm_coeff["K+"] = solve(zproc,Amin,Amax);
    std::cerr << "K+ Permeability Shift = " << perm_coeff["K+"] << std::endl;
}

void Cell:: RescaleEffectors()
{
    std::cerr << "Rescaling Effectors @Em=" << 1000 * Em << " mV" << std::endl;
    const double    zeta  = __Faraday__ * Em / ( __R__ * Temperature );
    
    //--------------------------------------------------------------------------
    // prepare solutions
    //--------------------------------------------------------------------------
    const solution &S_in   = inside(0.0);
    const solution &S_out  = outside(0.0);
    solution       &lambda = acc;
    solution       &rho    = tmp;
    //--------------------------------------------------------------------------
    // prepare leaks
    //--------------------------------------------------------------------------
    lambda.ldz();
    Leak(lambda, zeta, S_in, S_out);
    
    //--------------------------------------------------------------------------
    // prepare NaK
    //--------------------------------------------------------------------------
    {
        Effector    &NaK      = effectors["NaK"];
        rho.ldz();
        NaK.call(L, rho, 0.0, zeta, S_in);
        const double lambda_K = lambda["K+"];
        std::cerr << "\tlambda_K  = " << lambda_K << std::endl;
        if( lambda_K >= 0 )
            throw exception("Cell: invalid potassium concentrations");
        const double rho_NaK  = rho["K+"];
        std::cerr << "\trho_NaK   = " << rho_NaK << std::endl;
        NaK.factor = -lambda_K/rho_NaK;
        std::cerr << "\tx_NaK     = " << NaK.factor << std::endl;
    }
    
    //--------------------------------------------------------------------------
    // prepare AE2
    //--------------------------------------------------------------------------
    Effector &AE2 = effectors["AE2"];
    rho.ldz();
    AE2.call(L,rho,0.0,zeta,S_in);
    {
        const double lambda_Cl = lambda["Cl-"];
        std::cerr << "\tlambda_Cl = " << lambda_Cl << std::endl;
        if( lambda_Cl >= 0 )
            throw exception("Cell: invalid chloride concentrations");
        const double rho_AE = rho["Cl-"];
        std::cerr << "\trho_AE    = " << rho_AE << std::endl;
        AE2.factor = -lambda_Cl / rho_AE;
        std::cerr << "\tx_AE2     = " << AE2.factor << std::endl;
    }
    
    //--------------------------------------------------------------------------
    // prepare NHE
    //--------------------------------------------------------------------------
    Effector &NHE = effectors["NHE"];
    {
        rho.ldz();
        AE2.call(L,rho,0.0,zeta,S_in);
        const double rho_AE = rho["Cl-"];
        
        rho.ldz();
        NHE.call(L,rho,0.0,zeta,S_in);
        const double rho_NHE = rho["Na+"];
        std::cerr << "\trho_NHE   = " << rho_AE << std::endl;
        NHE.factor = rho_AE/rho_NHE;
        std::cerr << "\tx_NHE     = " << NHE.factor << std::endl;
    }
    
    std::cerr << "\tratio = " << AE2.factor / NHE.factor << std::endl;
}


void Cell:: RescaleAllRates( double coeff ) throw()
{
    
    std::cerr << "Rescaling All Rates by " << coeff << std::endl;
    const library &lib = Library;
    //! rescale permeabilities
    for( species::iterator p = lib.begin(); p != lib.end(); ++p )
    {
        Permeability &Perm = (**p).get<Permeability>();
        Perm.factor *= coeff;
    }
    
    //! rescale effectors
    for( Effectors::iterator q = effectors.begin(); q != effectors.end(); ++q )
    {
        Effector &eff = **q;
        eff.factor *= coeff;
    }
}


